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Abstract. Detection of the number of signals impinging on an array of sensors is an 
important problem in signal and array processing. Most of the papers consider asymptotic 
expansions of the sample size n whereas the dimension p of the observations is kept small. 
In this paper, we consider the case of high dimension, when p is large compared to n, using 
I recent results of random matrix theory. We extend our results obtained in the paper |13| 

T-H to the case of equal signals, and compare our algorithm to the method of Kritchman & 

O Nadler M, M- 

O 
H 

B 



1. Introduction 

In the area of signal processing, detection of the number of sources impinging on an array 
of sensors in presence of noise is a known and well- investigated problem [5], |6], ^0\, [T6] . 
This detection is generally a first step preliminary to any further study such as estimation 
of parameters. 



^ The underlying statistical model also appears in other scientific fields. In economics and 

T^ij- psychological literature, this model is called a factor model where the number of factors 

^ (signals) has a primary importance [1], [15]. Similar models can be found in physics of 

^ mixture [H], [H] or population genetics. More recently and in a slightly more general 

^ set-up, the model is introduced as a spiked population model 

I Many methods for determining the number of signals have been developed, mostly based 

'T^ on the minimum description length (MDL), Bayesian model selection or Bayesian Informa- 

tion Criteria (BIG) p5]. Nevertheless, these methods are based on asymptotic expansions 
^ for large sample size and may not perform well when the dimension of the data p is large 

^ compared to the sample size n. To avoid this problem of high dimension, several meth- 

ods have been recently proposed using the random matrix theory, such as Harding [7] or 
Onatski [12] in economics, and Kritchman & Nadler in chemometrics literature [9] and ar- 
ray processing [10]. In [13], we have also introduced a new method based on recent results 
of [2] and [H] in random matrix theory. 
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In all the cited references, signals are assumed to have distinct strengths. However, we 
observe that when some of these signals strengths become close, the detection problem 
is more difficult and most of existing algorithms need to be modified. We refer this new 
situation as the case with possibly equal signals and its precise formulation is given in 
Section 3.2. The aim of this work is to extend our method [13] to this new case and to 
compare it with the method of Kritchman & Nadler which, surprisingly enough, still 
work here even though it was introduced initially for the standard situation of distinct 
signals. 

The paper is organized as follows. In Section 2, we introduce the model. In Section 3, 
we define the detection problem of possibly equal signals and present our solution. We 
establish its asymptotic consistency. In Section 4, we recall the algorithm of Kritchman & 
Nadler [TD]. Next we conduct simulations experiments to compare these two methods. 

2. Problem Formulation 

Assuming an array of p sensors we consider the following standard model 

go 

(1) x(t) = ^afcSfc(t) + o-n(t) 

k=l 

(2) = As{t) + an{t), 
where 



• s{t) = (si(t), . . . , Sgp(t))* G M^" are go random signals assumed to have zero mean, unit 
variance and be mutually uncorrelated; 

• A = (ai, . . . , Ogg) is the p x qo steering matrix of go linearly independent p-dimensional 

vectors; 

• cr G M is the unknown noise level, n{t) ~ A/'(0, Ip) is a p x 1 vector of additive noise, 
independent of s{t). 

In this case, the population covariance matrix S = cov(x(t)) of x{t) takes the diagonal 
form 

W*SW = a\ + diag(ai, . . . , a,,, 0, . . . , 0) 

where W is an unknown basis of MP and ai > a2 > ■ ■ ■ > > 0. The sample covariance 
matrix of the n p-dimensional i.i.d. signal vectors received at each time t, (xj = x(tj))i<j<„ 
is 




i=l 



Denote by A„,i > Xn,2 > ■ ■ ■ > its eigenvalues. Our aim is to estimate go on the basis 
of S„. For the moment, we assume that the noise level o"^ is known. If this is not the case, 
we will give a method in section 4.2 to estimate it. 
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3. Detection of the number of signals 

In this section, we first recall our previous result of [13] in the case of different signals. 
Next, we propose an extension of the detection algorithm to the case with possibly equal 
signals. The consistency of the extended algorithm is established. 

3.1. Previous work: detection in the case of different signals. We consider the case 
where the {ai)i<i<qo are all different, so there is go distinct signals. According to [13], the 
population covariance matrix E has the spectral representation 

W*EW = aMiagK,..., «;,!,...,!), 

with the a^'s equals to 

It is assumed in the sequel that p and n are related so that when n ^ +oo, p/n ^ c > 0. 
This hypothesis allows the case where p is large compared to the sample size n (high- 
dimensional case). 

Moreover, we assumed that a[ > ■ ■ ■ > a'^^^ > 1 + y/c, i.e all the signal strength a are 

greater than a'^^fc. For a 7^ 1, we define the function 

, , , ca 

(p\a) = a^ -. 

a — 1 

Baik and Silverstein [3] proved that, under a moment condition on x, for each k G 
{1, . . . , go} and almost surely. 

They also proved that for all 1 < i < L with a prefixed range L and almost surely. 

The estimation method of the number of signals go in [13J is based on a close inspection of 
the following differences between consecutive eigenvalues 

^n,j = A„,j — A„j+i, J > 1. 

Indeed, from the results quoted above it is easy to see that a.s. if j > go, 5„j — )■ 
while when j < go, (5nj tends to a positive limit. Thus it is possible to detect go from 
index-numbers j where 5„j become small. More precisely, the estimator is 

(3) gn = min{j G {1, . . . , s} : 5„j+i < ci„}, 

where s > go is a fixed number big enough, and dn is a threshold to be defined. In practice, 
the integer s should be thought as a preliminary bound on the number of possible signals. 
In [13], we proved the consistency of g„ providing that the threshold satisfies dn — )■ 0, 



r^l'^dn — )■ -|-oo, under the following assumption on the entries of x. 



4 



DAMIEN PASSEMIER AND JIAN-FENG YAO 



Assumption 1. The entries x* of tlie random vector x liave a symmetric law and a sub- 
exponential decay, that is there exists positive constants C, C such that, for all t > C, 

P(|x'| > t^') < e-\ 

3.2. Detection in the case with possibly equal signals. As said in Introduction, 
when some of signals have close strengths, detection algorithms need to be modified. More 
precisely, we adopt the following theoretic model where we have K different signal strengths 
aK and for each signal strength a^, we have signals to be detected, that is 

spec(S) = (ai,...,ai,a2,...,a2,---,a£r-'_j_ax, O,--- ,0) + 0-^(1,...,!) 

ni 712 "A- P-90 P 

with ni + ■ ■ ■ + nK = Qo- When the signal strengths are different, the difference between the 
corresponding eigenvalues of the sample covariance matrix will tends to a positive constant, 
whereas with two signals of equal strength this difference will tends to zero. This fact 
creates an ambiguity with those differences corresponding to the noise eigenvalues which 
also tend to zero. Nevertheless, the convergence of the 6n,i, for i > qo (noise) is faster (in 
Op(n~^/^)) than that of the 5n/s corresponding to two identical signals (in Op{n~^^'^)) as 
a consequence of Theorem 3.1 of Bai & Yao [2]. This allows us to use the same estimator 
(|3|, provided we use a new threshold dn- The precise asymptotic consistency is as follows. 

Theorem 1. Let (xj)(i<j<„) be n copies i.i.d. ofx which follows the model and satisfies 
Assumption [1} We suppose that the population covariance matrix S has K non null and 
non unit eigenvalues ai > ... > ax > cr'^y/c with respective multiplicity {nk)i<k<K 
(ni + ■ ■ ■ + riK = qo), and p — qo unit eigenvalues. Assume that ^ — ). c > when n — )■ +oo. 
Let {dn)n>o be a real sequence such that dn = o(n~^/^) and n^^'^dn +oo. Then the 
estimator qn is strongly consistent, i.e g„, — )■ go almost surely when n — )■ +oo. 

Notice that the only modification of our estimator comparing to the different signals 
case is a new condition dn = o(?t,~^/^) on the convergence rate of dn- The proof of Theorem 
1 is postponed to Appendix. 

4. Method of Kritchman & Nadler with a noise estimation 

4.1. Algorithm of Kritchman & Nadler. In their paper [9] and [10], these authors 
develop a method based on another theorem from random matrix theory to detect the 
number of signals. 

In the absence of signals, nS„ follows a Wishart distribution with parameters n,p. In 
this case, Johnstone P] gave the asymptotic distribution of the largest eigenvalue of S„. 
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Proposition 1. Let S„ be the sample covariance matrix ofn vectors distributed as A/'(0, a^lp), 
and Xni > A„2 > ■■■ > ^np be its eigenvalues. Then, when n — )■ +00, such that 
2^00 

n 

where b = (l + i/c)^, (3n,p = (l + a/^) ^1 + y'^^ '^^'^ -^1 ^■^ Tracy-Widom distribution 
of order 1. 

Assuming the variance is known. To distinguish a signal eigenvalue A from a noise 
one at an asymptotic significance level 7, their idea is to check whether 

(4) Kk><yH^<i) + b 



n 

where 5(7) verifies ^1(5(7)) = 1 — 7 and can be found by inverting the Tracy-Widom 
distribution. This distribution has no explicit expression, but can be computed from a 
solution of a second order Painleve ordinary differential equation. Their estimator is based 
on a sequence of nested hypothesis tests of the following form: for A; = 1, 2, ... , min(p, n) — 1, 

"Hq: qo > k vs. l-Li. < k — \ . 

For each value of fc, they test the likelihood of the /c-th eigenvalue A„,fc as arising from a 
signal or from noise as ([i]). If Q is satisfied, I-Lq is accepted and k is increased by one. The 
procedure stops once an instance of 1-Lo is rejected and the number of signals is estimated 
to be g„ = /c — 1. Formally, their estimator is defined by 

qn = argmin ^A„,fc < ^ ^^^^^"^ ^(7) + ^ ~ ^• 

We refer this as the KN estimator. The authors proved the strong consistency of their 
algorithm as n — i- +00 with fixed p, by replacing the fixed confidence level 7 with a 
sample-size dependent one 7„, where 7^ — )■ sufficiently slow as n — > +00. They also 
proved that limp,„^+oo P {qn > qo) = 1- 

4.2. Estimation of the noise level. When the noise level o"^ is unknown, an estimation 
is needed. In [iSl, we used an algorithm based on the maximum likelihood estimate 

1 

p — qo ■ , 

As it is explained in [9] and [10], this estimator has a negative bias. Hence the authors 
developed an improved estimator with a smaller bias. We will use this improved estimator 
of noise level in our simulations for both estimator g„ and g„. 
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5. Simulation experiments and comparison 

In this section we compare by simulation our signal detector (PY) to the Kritchman & 
Nadler's one (KN). In their papers [S] and [TU], the authors compare the estimator KN with 
some other standard estimators in the signal processing literature, based on the minimum 
description length (MDL), Bayesian information criterion (BIG) and Akaike information 
criterion (AIC) [16j. In mostly studied cases, the estimator KN performs better. Thus we 
decided to consider only this detector for comparison. 

After several experiments, we found that the following version of our detector signifi- 
cantly improves the detection performance. Indeed, instead of stop once one difference 5k 
is below the threshold dn (see ([3])), the modified estimator now stops when two consecutive 
differences 5k and 5k+i are both below dn- More precisely, we set 

(5) g* = min{j G {1, . . . , s} : < dn and 5n,j+2 < dn]. 

It is easy to see that the proof for the consistency of g„ applies equally to g* under the 
same conditions as in Theorem 1. 

It remains to choose a threshold sequence dn to be used for our estimator g*. As argued 
in [T3], we use a sequence dn of the form Cn~'^^^y/2 log log n, where C is a constant to be 
adjusted for each case. In all simulations, we consider 500 independent replications and 
take = 1. We give a value of 7 = 0.5% to the false alarm rate of the estimator KN, as 
suggested in [10] and use their algorithm available at the author's homepage. 

Table 1 gives a summary of parameters in our simulation experiments. There are two 
sets of experiments for the comparison of the KN estimator with the ours g* (PY). In 
the first one (Figures 1-4 in the Table 1), signals have different strengths and these ex- 
periments extend and complete few results already reported in [13]. The second set of 
experiments (Figures 5-8 in Table 1) addresses the new situation where some signals have 
equal strengths. The last experiment (Figure 9) considers the case of no signal. 

5.1. Case of different signals. In Figure 1, we consider the case of a single signal with 
strength a, and we analyze the probability of misdetection as a function of the signal a, 
for (p,n) = (200,800), c = 0.25 and (p,n) = (2000,500), c = 4. We set C = 5.5 for the 
first case and C = 9 for the second case. The noise level cr^ = 1 is given to both estimators. 
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Table 1 . Summary of parameters used in the simulation experiments. 



Figure 


Signal 


Model 


Signals 




Fixed parameters 




Varying 


number 


strengths 


number 


strengths 


{P, n) 


c 


a' = l 


C 


parameters 


1 


Different 




(a) 


(200,800) 
(2000, 500) 


0.25 

4 


Given 


5.5 

9 


a 


2 


Different 


A 
B 


(6,5) 
(10,5) 




10 


Given 


11 


n 


3 


Different 


B 


(6,5) 




10 


To be estimated 


11 


n 


4 


Different 


C 
D 


(1.5) 
(2.5,1.5) 




1 


Given 


5 


n 


5 


Possibly 
equal 


E 
F 


(a, a, 5) 
(a, a, 15) 


(200,800) 
(2000, 500) 


0.25 

4 


Given 


6 

9.9 


a 


6 


Possibly 
equal 


G 
H 


(6,5,5) 
(10,5,5) 




10 


Given 


9.9 


n 


7 


Pes. equal 


H 


(10,5,5) 




10 


To be estimated 


9.9 


n 


8 


Possibly 
equal 


I 

J 


(1.5,1.5) 
(2.5,1.5,1.5) 




1 


Given 


4 


n 


9 


No signals 


K 


No signals 




1 

10 


Given 


8 

15 


n 




The two estimators have similar performance: the levels of detection are the same, and 
fit with the theory {^/c = 0.5 for the first case, and 2 for the second). 

In Figure 2, we consider two models with two signals (go = 2): 

• Model A: (ai,a2) = (6,5); 

• Model B: (ai,a2) = (10,5). 
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The detection is harder in Model A as the signal strengths are closer. We fix c = 10 
{p ^ n), and we plot the misdetection rates against the sample size n. Here C = 11. 
Again, o"^ = 1 is given to both estimators. 




As in Figure 1, the performances of the two estimator are close. However the estimator 
PY is slightly better for moderate values of n (n < 400) while the estimator KN has a 
slightly better performance for larger n. 

In Figure 3, we keep the same settings as in the previous simulation for Model B but 
with an unknown noise level cx^ to be estimated either. 




n 



Figure 3. Misdetection rates as a function of n for a ~ (10, 5) with estimated. 



Compared to Figure 2 (Model B), the estimation of o"^ does not affect the two estimators 
significantly. Both estimators seem robust against the unknown noise level. Note that we 
again observe a different hierarchy before and after n ~ 400. 
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Figure 4 considers two cases with c = 1 and a given noise level cr^ = 1: 

Model C: (a) = (1.5); 

Model D: (ai.as) = (2.5,1.5). 



Model C 



Model D 




— I — 

200 




— I — 

200 



Figure 4. Misdetection rates as a function of n for a = (1.5) (left) and a = (2.5, 1.5) (right). 

This experiment is designed with signal strengths getting closer to the critical value 
a/c = 1. This detection becomes more difficult and as expected both methods will have 
higher misdetection rates. Here we used C = 5. Meanwhile as displayed in Figure 4, 
our algorithm have a lower misdetection rate in almost all cases in both models, with an 
improvement ranging from 10% to 30% for moderate sample sizes n < 400. 



5.2. Case with equal signals. We keep the same parameters as in the previous section 
and only change the signal strengths. In Figure 5, we consider 

• Model E: (01,02,03) = (a, a, 5), < a < 2.5; 

• Model F: (ai, a2, as) = {a, a, 15), < a < 8. 

with (p, n) = (200, 800) for the Model E and {p, n) = (2000, 500) for the Model F. Here 
go = 3, C = 6 for Model E and C = 9.9 for Model F. 

This figure is to be compared to Figure 1 for different signals. Adding multiplicity only 
slightly increases the level of detection. As previously, the two estimators have similar 
performance. 
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Model E, (p=200,n=800) Model F, (p=2000,n=500) 




0.0 



—1 — 

0.5 



—\ — 

1.0 



— I— 

1.5 









































-0- KN 
-B- PY 





































2.0 



2.5 



Figure 5. Misdetection rates as a function of signal strength for {p,n) = (200,800) 
(left) and {p,n) = (2000,500) (right). 

In Figure 6, we consider two models analog to Model A and B with three signals (go = 3): 

Model G: {ai, 02,02) = (6,5,5); 
Model H: (01,02,02) = (10,5,5). 



Model G 



Model H 
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Figure 6. Misdetection rates as a function of n for a 

a = (10,5,5) (right). 



3, 5, 5) (left) and 



Again we fix c = 10 and we plot the probability of misdetection against the sample size 
n. Here C = 9.9 and cr^ is given. Comparing to the different signal strengths case (Figure 
2), the two estimators have significantly higher error rates. Nevertheless, the estimator 
PY shows superior detection performance for n < 500 (up to 20% less error): adding an 
equal signal affect more the performance of the estimator KN, but both estimators remain 
asymptotically consistent. If we compare Model G and Model H, a smaller spacing between 
the two first signals gives only a slightly degradation of correct detection. 
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In Figure 7, we keep the same settings as in the previous simulation for Model H but 
with an unknown noise level to be estimated either. 



Model H, estimated 




— I 1 1 1 r 

200 400 600 800 1000 



Figure 7. Misdetection rates as a function of n for a = (10, 5, 5) with cr^ estimated. 

Compared to Figure 6 (Model H), the estimation of cr^ does not affect the two estimators 
significantly. Both estimators seem robust against the unknown noise level. Note that we 
again observe a different hierarchy before and after n ~ 500. 

Figure 8 considers two cases with c = 1, and again = 1 is given: 

• Model I: (a, a) = (1.5,1.5); 

• Model J: (ai, aa, ^2) = (2.5, 1.5, 1.5). 

Model I Model J 




Figure 8. Misdetection rates as a function of n for a = (1.5, 1.5) (left) and 

a = (2.5,1.5,1.5) (right). 

Here we used C = 4. As explained in the previous section, this situation is more difficult 
and this causes a degradation in detection performance. The difference between the two 
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algorithms is higher than in the previous cases: the estimator PY performs better, up to 
40%. However, the convergence of both algorithms is quite slow. 

5.3. Case of no signal. In Figure 9 we examine the performance of the two algorithms 
in the case of no signal (Model K). The cases of c = 1 and c = 10 with cj^ = 1 given are 
considered. 



Model K, c=1 



Model K, 0=10 
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400 



— I — 
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— I — 
200 



— I — 

400 



— I — 

600 



Figure 9. Misdetection rates as a function of n in the case of no signals for c — 1 

(left) and c = 10 (right). 



We chose C = 8 for the first case and C = 15 for the second case. In both situations, 
the misdetection rates of the two estimators are similar and low (less than 1%). 



6. Concluding remarks 



In this paper we have considered the signal detection problem in the high-dimensional 
case. When some signals have close or even equal strengths, the detection becomes harder 
and existing algorithm need to be re-examined or corrected. In this spirit, we have proposed 
a new version of our previous algorithm. Its asymptotic consistency is established. It 
becomes unavoidable to compare our algorithm to an existing competitor proposed by 
Kritchman & Nadler (KN, |9], [in])- From our extensive simulation experiments in various 
scenarios, we observe that overall our detector has smaller misdetection rates, especially in 
cases with close and relatively low signal strengths (Figure 4 and 8) or more generally for 
almost all the cases provided that the sample size n is moderately large {n < 400 or 500). 

Another important conclusion concerns the influence of the tuning parameter, C for our 
algorithm and the false alarm rate 7 for KN. The main drawback of our algorithm is its lack 
of robustness with regard to the value of C. The experiments reported here are obtained 
with a finely-turned value of C and this value varies from case to case. How to overcome 
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this drawback, e.g by a data-adapted C remains an open problem. By comparison, the 
KN detector is remarkably robust and a single value of 7 = 0.5% was used in all the 
experiments. 

Appendix 

In the sequel, we will assume that o"^ = 1 (If it is not the case, we consider ^^). For 
the proof, we need two theorems. The first. Proposition [2| is a result of Bai and Yao [2] 
which derives a CLT for the n^-packed eigenvalues 

Vn[Xn,j - (Pia'k)], j G Jk 

where Jk = {s/t-i + !,•••, Sk}, Si = rii -\ h for 1 < z < i^. 

Proposition 2. Assume that the entries ofx satis/y E(||x*||^) < +00, a'^ > \^ ^ jor 
all 1 < j < K and have multiplicity ni, . . . ,nx respectively. Then as p, n ^ +00 so that 
- — )■ c, the nu- dimensional real vector 

Vn{Xn,j - (/'(a'fc); j e Jk} 

converges weakly to the distribution of the rik eigenvalues of a Gaussian random matrix 
whose covariance depend of a'f, and c. 

The second Proposition |3] is issued from the Proposition 5.8 of [1]: 

Proposition 3. Assume that the entries x* ofx have a symmetric law and a suh- exponential 
decay, that is there exists positive constants C, C such that, for allt > C, P(|x*| > t'") < e 
Then, for all 1 < i < L with a prefixed range L, 

2 

72.3 

where /3 = (1 + ^/c){l + \/(F^)^ . 

We also need the following lemma: 

Lemma 1. Let (X„)„>o be a sequence of positive random variables which converges weakly. 
Then for all real sequence (m„)„>o which converges to 0, 

F(X„ < Un) ^ 0. 

Proof. As (X„)„>o converges weakly, it exists a function G such that, for all v > 0, 
P(X„ < f ) — )■ G{v). Furthermore, cis y 0, it exists N eN such that for all n > N, 
Un < V. So P(X„ < Un) < P(X„ < v), and lim P(X„ < m„) < lim P(X„ < v) = 

n— >+oo n— >+oo 

G{v). Now we can take t> — )■ 0: as (X„)„>o is positive, G{v) 0. Consequently, 
P(Xn < Un) ^0. □ 



14 DAMIEN PASSEMIER AND JIAN-FENG YAO 

Proof, of Theorem [l] The proof is essentially the same as Theorem 3.1 in [13], except when 
the spikes are equal. We have 

{qn = go} = {go = min{j : 6n,j+i < dn}} 

= {Vj G {1, . . . , go}, Sn,j > dn} n {Sn,qo+l < dn} ■ 

Therefore 

P(gn = go) = P( fl {Sn,j > dn} n {5n,go+l < dn}] 

\i<j<go / 

= 1-P I U {Snj < dn} {6n,go+l > dn}] 

\i<i<go / 

<?o 

> H^n,j < dn) - P(5n,,o + l > dn) ■ 

Case of j = qq + 1. In this case, ^n.qo+i — -^n.go+i ~ ^n,go+2 (noise eigenvalues). As (i„ — )■ 
such that, n'^^^dn — )■ +oo, and by using Proposition [3] in the same manner as in the proof 
of Theorem 3.1 in [13], we have 

F(5„,go+i > dn) -> 0. 



Case of 1 < j < go- These indices correspond to the signal eigenvalues. 

• Let /i = {1 < / < go|card(J/) = 1} (simple signal) and I2 = — 1|^ G h}- For all 
j G Ji U /2, Sn,j corresponds to a consecutive difference of A„j issued from two different 
signals, so we can still use Proposition [2] and the proof of Theorem 3.1 in [T3] to show 
that 

F{6nj < dn) 0, Vj G Ji. 

• Let Jg = {1 < / < go - ^ ih U h)}- For all j G I3, it exists k e {1,. . . ,K} such 
that j G Jfc. By Proposition [2| X„ = \/n6n,j converges weakly. So by using Lemma [l] 
and that dn = o(n~^/^), we have 

P {5nJ <dn)=F {^/^5n,j < V^dn) ^ 0. 

• The case of j = go + 1 is considered as in [T3] . 

Conclusion. F{6n,qo+i > dn) — )• and Ylf=i ^i^^nj < dn) 0, therefore 
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P(?n = go) 1- 

n—^+oo 

□ 
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